module debug_utils
  implicit none

  contains

  !> Calculate the L2 norm of an array and print the result when rank is 0
  !>
  !> @param classname The class identifier (e.g., "Cdebug" or "Fdebug")
  !> @param TAG The identifier to print with the norm
  !> @param array The array to calculate the norm for
  !> @param size The size of the array
  !> @param rank The MPI rank (only rank 0 will print)
  !> @param step Optional step counter to make each TAG unique (default: 0)
  subroutine debug_print_norm(classname, TAG, array, size, rank, step)
    character(len=*), intent(in) :: classname
    character(len=*), intent(in) :: TAG
    real(kind=4), intent(in) :: array(*)
    integer, intent(in) :: size
    integer, intent(in) :: rank
    integer, intent(in), optional :: step
    
    real :: sum_squares
    real :: norm
    integer :: i
    integer :: step_val
    
    ! Use default step value of 0 if not provided
    if (present(step)) then
      step_val = step
    else
      step_val = 0
    end if
    
    if (rank == 0) then
      sum_squares = 0.0
      do i = 1, size
        sum_squares = sum_squares + array(i) * array(i)
      end do
      norm = sqrt(sum_squares)
      ! Use scientific notation with 5 decimal places for consistent output
      write(*,'(A,A,A,A,A,I0,A,ES12.5)') '[', trim(classname), '] ', trim(TAG), ' step=', step_val, ' norm=', norm
    end if
  end subroutine debug_print_norm
  
  !> Calculate the L2 norm of a double-precision array and print the result when rank is 0
  !>
  !> @param classname The class identifier (e.g., "Cdebug" or "Fdebug")
  !> @param TAG The identifier to print with the norm
  !> @param array The array to calculate the norm for
  !> @param size The size of the array
  !> @param rank The MPI rank (only rank 0 will print)
  !> @param step Optional step counter to make each TAG unique (default: 0)
  subroutine debug_print_norm_d(classname, TAG, array, size, rank, step)
    character(len=*), intent(in) :: classname
    character(len=*), intent(in) :: TAG
    real(kind=8), intent(in) :: array(*)
    integer, intent(in) :: size
    integer, intent(in) :: rank
    integer, intent(in), optional :: step
    
    real(kind=8) :: sum_squares
    real(kind=8) :: norm
    integer :: i
    integer :: step_val
    
    ! Use default step value of 0 if not provided
    if (present(step)) then
      step_val = step
    else
      step_val = 0
    end if
    
    if (rank == 0) then
      sum_squares = 0.0d0
      do i = 1, size
        sum_squares = sum_squares + array(i) * array(i)
      end do
      norm = sqrt(sum_squares)
      ! Use scientific notation with 12 decimal places for higher precision output
      write(*,'(A,A,A,A,A,I0,A,ES19.12)') '[', trim(classname), '] ', trim(TAG), ' step=', step_val, ' norm=', norm
    end if
  end subroutine debug_print_norm_d
  
end module debug_utils 